# updated 1/19/2015 to work with R 3.1.2 and lme4 1.1-7
# also commented out graphics and summaries
# ======================================================
# Now estimate boltival from July 2007 to July 2008
# ======================================================

bolt=read.csv("July 2007 census and biomass2 edited.csv",header = T,na.strings=c("NA","na"))
#converting variables in "bolt" data to numeric/categorical variables
bolt$comp = as.factor(bolt$comp)
bolt$herb = as.factor(bolt$herb)
bolt$species = as.factor(bolt$species)
bolt$bolt08 = as.factor(bolt$bolt)
# figure out which rows have the ugly "" factor levels
bolt[bolt$Comp=="",5:10] = NA
# relevel competition factor
bolt$Comp = factor(bolt$Comp,levels=c("Low","Medium","High"))
bolt$Species = factor(bolt$Species,levels=c("BT","TT"))
bolt$Herb = factor(bolt$Herb,levels=c("Ambient","Reduced"))
bolt = bolt[complete.cases(bolt[,5:10]),]

with(bolt,table(Species,survived08))

# Check for complete cases of the variables we need 
# no.leaves07, survived08, bolt08
# alot of missing values for no.leaves07 in bull thistle
with(bolt,table(Species,is.na(no.leaves07)))
bolt = bolt[complete.cases(bolt[,c("no.leaves07","survived08","bolt08")]),]
with(bolt,table(bolt08,survived08,Species)) 

# for compatibility, split into tall and bull datasets, discard plants that didn't survive
bolt.tall = bolt[bolt$Species=="TT" & bolt$survived08>0,]
bolt.bull = bolt[bolt$Species=="BT" & bolt$survived08>0,]
bolt.tall$plot = paste(bolt.tall$block,bolt.tall$plot,sep=".")
bolt.bull$plot = paste(bolt.bull$block,bolt.bull$plot,sep=".")

library(lme4)
###########################################################################################
## Fit boltival models for tall thistle
###########################################################################################
# now predict veg.bio07 for this dataset. 
# assumes that growth models have been fitted, check for appropriate models
if (!exists("tb.fits")) stop("fit growth models first!")

# now using model averaged predictions
nd = bolt.tall
nd$no.leaves08  = nd$no.leaves07
nd$bolt08 = factor(0,levels=levels(grow.tall$bolt08))
results = matrix(NA,nrow=nrow(bolt.tall),ncol=length(tb.fits))
for (i in 1:length(tb.fits)){
	# generate predictions from a lmer model using the data from inside the model, mostly. 
	fm = tb.fits[[i]]
  # only use block random effects, because some plots don't have growth data
	results[,i] = predict(fm,newdata=nd,re.form=~(1|block))
}
# model averaged prediction for veg.bio07 parameter
# make sure we have the right weights vector!!
bic = unlist(sapply(tb.fits,BIC)		)
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
bolt.tall[,"veg.bio07"] = apply(results,1,function(x)sum(x*wbic))

# fit global model with block and plot effects
blt.mm0 = glmer(bolt08~veg.bio07 + Comp + Herb + 
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|block) + (1|plot), family=binomial, data=bolt.tall)
# compare model with variation in herbivory among blocks
blt.mm1 = glmer(bolt08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(Herb|block) + (1|plot), family=binomial, data=bolt.tall)
blt.mm2 = glmer(bolt08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|block) + (Herb|plot), family=binomial, data=bolt.tall)

blt.mm3 = glmer(bolt08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|block) + (Comp|plot), family=binomial, data=bolt.tall)
blt.mm4 = glmer(bolt08~veg.bio07 + Comp + Herb + 
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(Comp|block) + (1|plot), family=binomial, data=bolt.tall)
# worse, try removing random effects
blt.mm5 = glmer(bolt08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|block), family=binomial, data=bolt.tall)
# NO! very bad, remove block?
blt.mm6 = glmer(bolt08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|plot), family=binomial, data=bolt.tall)
sapply(list(blt.mm0,blt.mm1,blt.mm2,blt.mm3,blt.mm4,blt.mm5,blt.mm6),BIC)
# OK, simple plot effect only

# OK try as BIC model selection to do multi-model inference
blt.models = list("~veg.bio07 + Herb + Comp + veg.bio07:Comp + veg.bio07:Herb + Herb:Comp",
		"~veg.bio07 + Herb + Comp + veg.bio07:Comp + veg.bio07:Herb",
		"~veg.bio07 + Herb + Comp + veg.bio07:Comp + Herb:Comp",
		"~veg.bio07 + Herb + Comp + veg.bio07:Herb + Herb:Comp",
		"~veg.bio07 + Herb + Comp + veg.bio07:Herb",
		"~veg.bio07 + Herb + Comp + Herb:Comp",
		"~veg.bio07 + Herb + Comp + veg.bio07:Comp",
		"~veg.bio07 + Herb + Comp",
		"~veg.bio07 + Herb +  veg.bio07:Herb",
		"~veg.bio07 + Herb",
		"~veg.bio07 + Comp + veg.bio07:Comp",
		"~veg.bio07 + Comp",
		"~Herb + Comp + Herb:Comp",
		"~Herb + Comp",
		"~Herb",
		"~Comp",
		"~veg.bio07",
		"~1")

blt.fits = lapply(blt.models,function(ff)glmer(paste("bolt08",ff,"+ (1|plot)"),data=bolt.tall,family=binomial))
bic = unlist(sapply(blt.fits,BIC))
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
data.frame(unlist(blt.models),dbic,wbic)[order(dbic),]
# model 12 has 88% of the weight, and model 17 has 4.5% - pretty much vegbio07 and comp
loglik=unlist(sapply(blt.fits,logLik))
blt.fits[[1]]
k=sapply(blt.fits,function(mm)length(fixef(mm)))+1
stat.res=data.frame(unlist(blt.models),k,loglik,dbic,wbic)[order(dbic),]
write.csv(stat.res,"bolting_res_TT.csv")

# make a plot
nd = data.frame(veg.bio07=rep(seq(-2.5,0.8,0.1),times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=68),levels=levels(bolt.tall$Comp)),
		Herb=factor(rep(rep(c("Ambient","Reduced"),3),each=34)))

results = matrix(NA,nrow=nrow(nd),ncol=length(blt.fits))
for (i in 1:length(blt.fits)){
	# generate predictions from a glmer model using the data from inside the model, mostly. 
	fm = blt.fits[[i]]
	results[,i] = predict(fm,newdata=nd,re.form=~0)
}
# model averaged prediction for veg.bio07 parameter
nd[,"bolt"] = apply(results,1,function(x)sum(x*wbic))

xyplot(bolt~veg.bio07|Comp+Herb,data=nd,type="l",
panel=function(...){
		panel.xyplot(...)
		pn = which.packet()
		#ltext(0,0.5,pn[1])
		pick = bolt.tall$Comp == levels(bolt.tall$Comp)[pn[1]] & 
			 bolt.tall$Herb == levels(bolt.tall$Herb)[pn[2]]
		ticks = bolt.tall[pick,c("veg.bio07","bolt08")]
		panel.rug(ticks[ticks$bolt08==0,1])
		panel.rug(ticks[ticks$bolt08==1,1],y=NULL, regular=FALSE)
	},
ylim=c(0,1))
#pull out the parameters of teh best model
fixef(blt.fits[[12]])


#------------------
coef.names = names(fixef(blt.fits[[1]]))
coef.results = matrix(NA,nrow=length(blt.fits),ncol=length(coef.names))
coef.se = matrix(NA,nrow=length(blt.fits),ncol=length(coef.names))
for (i in 1:length(blt.fits)){
	coef.results[i,] = sapply(coef.names,function(cf)fixef(blt.fits[[i]])[cf],USE.NAMES=FALSE)
	se = sqrt(diag(vcov(blt.fits[[i]])))
	names(se) = names(fixef(blt.fits[[i]]))
	coef.se[i,] = sapply(coef.names,function(cf)se[cf])
}

coef.avg = apply(coef.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(coef.results,2,coef.avg)
t2 = sqrt(coef.se^2 + t1^2)
coef.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

nd.intercept = data.frame(veg.bio07=rep(0,times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=2),levels=levels(bolt.tall$Comp)),
		Herb=factor(rep(c("Ambient","Reduced"),3)))
mm.intercept = model.matrix(as.formula(blt.models[[1]]),data=nd.intercept)

nd.slope = data.frame(veg.bio07=rep(1,times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=2),levels=levels(bolt.tall$Comp)),
		Herb=factor(rep(c("Ambient","Reduced"),3)))

mm.slope = model.matrix(as.formula(blt.models[[1]]),data=nd.slope)

mm.slope[,c(1,3:5,9:10)] = 0

int.results = array(NA,dim=c(length(blt.fits),nrow(mm.intercept)))
int.se = matrix(NA,nrow=length(blt.fits),ncol=nrow(mm.intercept))
slope.results = array(NA,dim=c(length(blt.fits),nrow(mm.slope)))
slope.se = matrix(NA,nrow=length(blt.fits),ncol=nrow(mm.intercept))

for (i in 1:length(blt.fits)){
	fef = rep(0,times=length(coef.names))
	names(fef) = coef.names
	vcf = matrix(0,nrow=length(coef.names),ncol=length(coef.names))
	rownames(vcf) = coef.names
	colnames(vcf) = coef.names

	fe = fixef(blt.fits[[i]])
	fef[names(fe)] = fe
	vcf[names(fe),names(fe)] = as.matrix(vcov(blt.fits[[i]]))
	int.results[i,] = mm.intercept %*% fef
	slope.results[i,] = mm.slope %*% fef
	for (j in 1:nrow(mm.intercept)){
		pick = which(mm.intercept[j,]>0)
		int.se[i,j] = sqrt(sum(vcf[pick,pick]))
		pick = which(mm.slope[j,]>0)
		slope.se[i,j] = sqrt(sum(vcf[pick,pick])) # make it the se not the variance
	}
	
}

int.avg = apply(int.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(int.results,2,int.avg)
t2 = sqrt(int.se^2 + t1^2)
int.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

slope.avg = apply(slope.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(slope.results,2,slope.avg)
t2 = sqrt(slope.se^2 + t1^2)
slope.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

for (i in 1:6) {IPM.parameters$bolt.intercept[IPM.parameters$species=="TT" & IPM.parameters$comp==nd.intercept[i,"Comp"] & IPM.parameters$herb==nd.intercept[i,"Herb"]]=int.avg[i]}
for (i in 1:6) {IPM.parameters$bolt.intercept.se[IPM.parameters$species=="TT" & IPM.parameters$comp==nd.intercept[i,"Comp"] & IPM.parameters$herb==nd.intercept[i,"Herb"]]=int.avg.se[i]}
for (i in 1:6) {IPM.parameters$bolt.slope[IPM.parameters$species=="TT" & IPM.parameters$comp==nd.slope[i,"Comp"] & IPM.parameters$herb==nd.intercept[i,"Herb"]]=slope.avg[i]}
for (i in 1:6) {IPM.parameters$bolt.slope.se[IPM.parameters$species=="TT" & IPM.parameters$comp==nd.slope[i,"Comp"] & IPM.parameters$herb==nd.intercept[i,"Herb"]]=slope.avg.se[i]}







###########################################################################################
## Fit boltival models for Bull thistle
###########################################################################################
# now predict veg.bio07 for this dataset. 
# assumes that growth models have been fitted, check for appropriate models
if (!exists("bb.fits")) stop("fit growth models first!")

# now using model averaged predictions
nd = bolt.bull
nd$no.leaves08  = nd$no.leaves07
nd$bolt08 = factor(0,levels=levels(grow.bull$bolt08))
results = matrix(NA,nrow=nrow(bolt.bull),ncol=length(bb.fits))
for (i in 1:length(bb.fits)){
	# generate predictions from a lmer model using the data from inside the model, mostly. 
	fm = bb.fits[[i]]
	results[,i] = predict(fm,newdata=nd,re.form=~(1|block))
}
# model averaged prediction for veg.bio07 parameter
# make sure we have the right weights
bic = unlist(sapply(bb.fits,BIC)		)
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
bolt.bull[,"veg.bio07"] = apply(results,1,function(x)sum(x*wbic))

# fit global model with block and plot effects
blb.mm0 = glmer(bolt08~veg.bio07 + Comp + Herb + 
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|block) + (1|plot), family=binomial, data=bolt.bull)
# compare model with variation in herbivory among blocks
blb.mm1 = glmer(bolt08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(Herb|block) + (1|plot), family=binomial, data=bolt.bull)
blb.mm2 = glmer(bolt08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|block) + (Herb|plot), family=binomial, data=bolt.bull)

blb.mm3 = glmer(bolt08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|block) + (Comp|plot), family=binomial, data=bolt.bull)
blb.mm4 = glmer(bolt08~veg.bio07 + Comp + Herb + 
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(Comp|block) + (1|plot), family=binomial, data=bolt.bull)
# worse, try removing random effects
blb.mm5 = glmer(bolt08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|block), family=binomial, data=bolt.bull)
# NO! very bad, remove block?
blb.mm6 = glmer(bolt08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|plot), family=binomial, data=bolt.bull)
sapply(list(blb.mm0,blb.mm1,blb.mm2,blb.mm3,blb.mm4,blb.mm5,blb.mm6),BIC)
# OK, simple plot effect only

# OK try as BIC model selection to do multi-model inference
blb.models = list("~veg.bio07 + Herb + Comp + veg.bio07:Comp + veg.bio07:Herb + Herb:Comp",
		"~veg.bio07 + Herb + Comp + veg.bio07:Comp + veg.bio07:Herb",
		"~veg.bio07 + Herb + Comp + veg.bio07:Comp + Herb:Comp",
		"~veg.bio07 + Herb + Comp + veg.bio07:Herb + Herb:Comp",
		"~veg.bio07 + Herb + Comp + veg.bio07:Herb",
		"~veg.bio07 + Herb + Comp + Herb:Comp",
		"~veg.bio07 + Herb + Comp + veg.bio07:Comp",
		"~veg.bio07 + Herb + Comp",
		"~veg.bio07 + Herb +  veg.bio07:Herb",
		"~veg.bio07 + Herb",
		"~veg.bio07 + Comp + veg.bio07:Comp",
		"~veg.bio07 + Comp",
		"~Herb + Comp + Herb:Comp",
		"~Herb + Comp",
		"~Herb",
		"~Comp",
		"~veg.bio07",
		"~1")

blb.fits = lapply(blb.models,function(ff)glmer(paste("bolt08",ff,"+ (1|plot)"),data=bolt.bull,family=binomial))
bic = unlist(sapply(blb.fits,BIC))
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
data.frame(unlist(blb.models),dbic,wbic)[order(dbic),]
# model 17 has 83% of the weight, and model 12 has 9%, then come models 10 and 9 
loglik=unlist(sapply(blb.fits,logLik))
k=sapply(blb.fits,function(mm)length(fixef(mm)))+1
stat.res=data.frame(unlist(blb.models),k,loglik,dbic,wbic)[order(dbic),]
write.csv(stat.res,"bolting_res_BT.csv")

# make a plot
nd = data.frame(veg.bio07=rep(seq(-2.5,0.8,0.1),times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=68),levels=levels(bolt.bull$Comp)),
		Herb=factor(rep(rep(c("Ambient","Reduced"),3),each=34)))

results = matrix(NA,nrow=nrow(nd),ncol=length(blb.fits))
for (i in 1:length(blb.fits)){
	# generate predictions from a glmer model using the data from inside the model, mostly. 
	fm = blb.fits[[i]]
	results[,i] = predict(fm,newdata=nd,re.form=~0)
}
# model averaged prediction for veg.bio07 parameter
nd[,"bolt"] = apply(results,1,function(x)sum(x*wbic))

xyplot(bolt~veg.bio07|Comp+Herb,data=nd,type="l",
panel=function(...){
		panel.xyplot(...)
		pn = which.packet()
		#ltext(0,0.5,pn[1])
		pick = bolt.bull$Comp == levels(bolt.bull$Comp)[pn[1]] & 
			 bolt.bull$Herb == levels(bolt.bull$Herb)[pn[2]]
		ticks = bolt.bull[pick,c("veg.bio07","bolt08")]
		panel.rug(ticks[ticks$bolt08==0,1])
		panel.rug(ticks[ticks$bolt08==1,1],y=NULL, regular=FALSE)
	},
ylim=c(0,1))
fixef(blb.fits[[17]])


#------------------
coef.names = names(fixef(blb.fits[[1]]))
coef.results = matrix(NA,nrow=length(blb.fits),ncol=length(coef.names))
coef.se = matrix(NA,nrow=length(blb.fits),ncol=length(coef.names))
for (i in 1:length(blb.fits)){
	coef.results[i,] = sapply(coef.names,function(cf)fixef(blb.fits[[i]])[cf],USE.NAMES=FALSE)
	se = sqrt(diag(vcov(blb.fits[[i]])))
	names(se) = names(fixef(blb.fits[[i]]))
	coef.se[i,] = sapply(coef.names,function(cf)se[cf])
}

coef.avg = apply(coef.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(coef.results,2,coef.avg)
t2 = sqrt(coef.se^2 + t1^2)
coef.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

nd.intercept = data.frame(veg.bio07=rep(0,times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=2),levels=levels(bolt.bull$Comp)),
		Herb=factor(rep(c("Ambient","Reduced"),3)))
mm.intercept = model.matrix(as.formula(blb.models[[1]]),data=nd.intercept)

nd.slope = data.frame(veg.bio07=rep(1,times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=2),levels=levels(bolt.bull$Comp)),
		Herb=factor(rep(c("Ambient","Reduced"),3)))

mm.slope = model.matrix(as.formula(blb.models[[1]]),data=nd.slope)

mm.slope[,c(1,3:5,9:10)] = 0

int.results = array(NA,dim=c(length(blb.fits),nrow(mm.intercept)))
int.se = matrix(NA,nrow=length(blb.fits),ncol=nrow(mm.intercept))
slope.results = array(NA,dim=c(length(blb.fits),nrow(mm.slope)))
slope.se = matrix(NA,nrow=length(blb.fits),ncol=nrow(mm.intercept))

for (i in 1:length(blb.fits)){
	fef = rep(0,times=length(coef.names))
	names(fef) = coef.names
	vcf = matrix(0,nrow=length(coef.names),ncol=length(coef.names))
	rownames(vcf) = coef.names
	colnames(vcf) = coef.names

	fe = fixef(blb.fits[[i]])
	fef[names(fe)] = fe
	vcf[names(fe),names(fe)] = as.matrix(vcov(blb.fits[[i]]))
	int.results[i,] = mm.intercept %*% fef
	slope.results[i,] = mm.slope %*% fef
	for (j in 1:nrow(mm.intercept)){
		pick = which(mm.intercept[j,]>0)
		int.se[i,j] = sqrt(sum(vcf[pick,pick]))
		pick = which(mm.slope[j,]>0)
		slope.se[i,j] = sqrt(sum(vcf[pick,pick])) # make it the se not the variance
	}
	
}

int.avg = apply(int.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(int.results,2,int.avg)
t2 = sqrt(int.se^2 + t1^2)
int.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

slope.avg = apply(slope.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(slope.results,2,slope.avg)
t2 = sqrt(slope.se^2 + t1^2)
slope.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))


for (i in 1:6) {IPM.parameters$bolt.intercept[IPM.parameters$species=="BT" & IPM.parameters$comp==nd.intercept[i,"Comp"] & IPM.parameters$herb==nd.intercept[i,"Herb"]]=int.avg[i]}
for (i in 1:6) {IPM.parameters$bolt.intercept.se[IPM.parameters$species=="BT" & IPM.parameters$comp==nd.intercept[i,"Comp"] & IPM.parameters$herb==nd.intercept[i,"Herb"]]=int.avg.se[i]}
for (i in 1:6) {IPM.parameters$bolt.slope[IPM.parameters$species=="BT" & IPM.parameters$comp==nd.slope[i,"Comp"] & IPM.parameters$herb==nd.intercept[i,"Herb"]]=slope.avg[i]}
for (i in 1:6) {IPM.parameters$bolt.slope.se[IPM.parameters$species=="BT" & IPM.parameters$comp==nd.slope[i,"Comp"] & IPM.parameters$herb==nd.intercept[i,"Herb"]]=slope.avg.se[i]}


write.csv(IPM.parameters,"IPM.Parameters.csv",row.names = FALSE)